/** Author: David Powell
This file produces Figure A.15
Overdoses are suppressed in data file
**/


clear all
set more off
set mat 800
set seed 8721

global dir "/jules/b/dpowell"
global DATA "${dir}/purdue/replication/DATA"
global OUTPUT "${dir}/purdue/replication/output"

use countydata

gen overdose_rate=100000*(overdose/totpop)

tab year, gen(ttt) nof
drop ttt13
foreach nnn of numlist 1983/1994 1996/2017 {
	gen tt_`nnn'=nontripl*(year==`nnn')
}

egen ggg1=group(stf county)

**balanced panel
drop if totpop==.
bys ggg1: gen numobs=_N
keep if numobs==35
drop numobs 



gen large=(ruralcode<4)
gen verylarge=(ruralcode<1)


gen beta1=0
gen low1=0
gen high1=0
xtivreg2 overdose_rate tt_* ttt* if large==1 [aw=totpop], cluster(stf)  i(ggg1) fe
foreach nnn2 of numlist 1983/1994 1996/2017 {
	replace beta1=_b[tt_`nnn2'] if year==`nnn2'
	boottest tt_`nnn2', boottype(wild) weighttype(webb)  reps(9999) seed(23857389)
	matrix A=r(CI)
	replace low1=A[1,1] if year==`nnn2'
	replace high1=A[1,2] if year==`nnn2'		
}


gen beta2=0
gen low2=0
gen high2=0
xtivreg2 overdose_rate tt_* ttt* if verylarge==1 [aw=totpop], cluster(stf)  i(ggg1) fe
foreach nnn2 of numlist 1983/1994 1996/2017 {
	replace beta2=_b[tt_`nnn2'] if year==`nnn2'
	boottest tt_`nnn2', boottype(wild) weighttype(webb)  reps(9999) seed(23857389)
	matrix A=r(CI)
	replace low2=A[1,1] if year==`nnn2'
	replace high2=A[1,2] if year==`nnn2'		
}

bys year: keep if _n==1
# delimit ;

twoway(rarea low1 high1 year,  color(gs14))  (connected beta1 year, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1)) ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figA15a.eps, replace;

twoway(rarea low2 high2 year,  color(gs14))  (connected beta2 year, color(black) msize(vtiny) sort yaxis(1)) 
, graphregion(color(white)) xline(1996) yline(0, lpattern(dash) lcolor(gs6)) xlabel(1983(2)2017, angle(45))  legend(label(2 "Estimate") label(1 "95% Confidence Interval")  rows(1) order(2 1)) ylabel(,nogrid)
scheme(s2mono) xtitle("") ytitle("Coefficient Estimate" ,  axis(1)) ;
gr export ${OUTPUT}/figA15b.eps, replace;

# delimit cr
